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Abstract 

Technical work during September-December 1990 consisted of 1) analyzing 
HST point source images obtained from JPL, 2) retrieving phase information 
from the images by a direct (noniterative) technique, and 3) characterizing the 
wavefront aberration due to errors in the Hubble Space Telescope (HST) 
mirrors, in a preliminary manner. This work was in support of JPL design of 
compensating optics for the next generation wide-field planetary camera on 
HST. This digital technique for phase retrieval from pairs of defocused images, 
is based on the energy transport equation between these image planes. In 
addition, an end-to-end wave optics routine, based on the JPL Code V 
prescription of the unaberrated HST and WFPC, was derived for output of the 
reference phase front when mirror error is absent. Also, the Roddier routine 
unwrapped the retrieved phase by inserting the required jumps of ±2n radians 
for the sake of smoothness. A least-squares fitting routine, insensitive to phase 
unwrapping, but nonlinear, was used to obtain estimates of the Zernike 
polynomial coefficients that describe the aberration. The phase results were 
close to, but higher than, the expected error in conic constant of the primary 
mirror suggested by the fossil evidence. The analysis of aberration contributed 
by the camera itself could be responsible for the small discrepancy, but was not 
verified by analysis. The wavefront Zernike coefficient Z u (in HST OTA 
Handbook convention) was found to be -0.27 and -0.28 |im for the HI -II 
(HARP 1 A data)and 01 -PI (HARP IB data) pairs of images, respectively. The 

standard deviation error in these quantities was approximately ±0.005, due to 
the iterations with the nonlinear fitting routine for Zernike coefficients. The 
implied conic constant for the primary mirror was -1 .01 41 and -1 .01 46 ±0.0002 
compared to the design value of -1.0022985. 

New Technology . No new technology was identified. 
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1.0 Introduction 


The Hubble Space Telescope (HST) has not performed as well as designed 
owing to an imperfect mirror (or mirrors). The Jet Propulsion Laboratory (JPL) is 
to define the mirror imperfections and to plan a hardware fix in their Wide Field 
Planetary Camera (WF/PC) to compensate for the aberrations. The work 
reported here supported JPL in defining the HST aberration by means of a 
direct phase retrieval technique. This technique operated on digital images of 
stars to obtain the shape of the optical wavefront in the exit pupil plane. This 
wavefront distortion was defined in classical terms such as spherical aberration, 
coma, astigmatism, and so on. This wavefront analysis also would lead to a 
Code V prescription, for example, of optical elements of the WF/PC that would 
be needed to compensate for the aberrations. 

Early on, it was recognized that spherical aberration was present in large 
amounts. This led to the suspicion of an incorrect conic constant, which 
specifies the hyperboloid shape of the primary mirror. The purpose of the 
phase retrieval work was to quantify this suspicion, and to determine if other 
errors existed in the primary mirror, other errors might be in the secondary 
mirror and the camera. A strategy of using multiple off-axis images for 
distinguishing secondary from primary mirror errors was voiced by various 
participants at the first HARP Image Inversion Workshop. To date little work 
toward this end has been done due to a lack of images at sufficiently large field 
angles. 

The imperfections in the camera might be another story. There is a difference in 
conic constant error estimates from the fossil evidence (interferograms and 
notebooks at Perkin-Elmer, now Hughes Danbury Optical Systems), and from 
the Zernike polynomial coefficient in the phase retrieval effort. This, albeit small, 
difference of conic constant might be attributed to the camera. End-to-end 
modeling that includes the camera, with its distinguishing obscurations, has 
been part of the phase retrieval effort. 

Much of the phase retrieval work by others consisted of 1) least-squares fitting 
of images generated by OTA-camera models with image data, 2) iterative 
automatic fitting by repeated single or double Fourier transforming between 
pairs or triplets of image planes and their corresponding pupil planes, and 3) 
neural network Zernike polynomial output after training with simulated images 
of known Zernike polynomial input. The work reported here was a direct 
(noniterative) phase retrieval technique that solves a nonlinear differential 
equation to complete the field information in the image plane so that a Fourier 
transform retrieves the pupil plane phase that was sought. 

According to PC data, the Airy disk of the point spread function (PSF) would 
cover a few pixels if the OTA operated as designed. The present condition 
seemed to result in a peak nearly as sharp but with significant intensity in rings 
and spokes around the peak. 
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2.0 Technical Purpose 

The work consisted of the following three technical tasks for defining the Hubble 
Space Telescope (HST) aberration by means of the direct (noniterative) phase 
retrieval technique developed at TRW. This technique was based on the 
Gonsalves 'phase diversity' algorithm, which processed pairs of images that 
differed by a small amount of defocus. 

2.1 Task 1 . Analysis of JPL Star Images 

Digital images were input to the phase retrieval code. The code was modified 
as follows: 1) to read-in the image data, 2) to define a limited field of view 
containing each star image, and 3) to register (align) the focus/defocus image 
pairs. 

Simulation of the image intensity distribution for these different cases of 
wavelength filter and secondary mirror positions was done for two purposes: 1 ) 
to allow tests of the phase retrieval algorithm for known amounts of aberration 
and 2) to enable interpolation between pixels of data for more accurate phase 
retrieval. 

Only cases of star objects lying on axis were analyzed. 

2.2 Task 2. Implementation of phase retrieval code 

Wave optics of the HST optical train was simulated, as needed for phase 
retrieval purposes. The code was modified 1) for the effects of apertures, 
obscurations, etc., and 2) for the mirror displacement aberrations. The phase 
retrieval code was executed on a VAX 8650 system and a 386-PC (25 MHz, 1 6 
Mbytes RAM) 1 ) to produce Zernike polynomial fits to the wavefront, and 2) to 
estimate error limits on the phase retrieval accuracy. 

2.3 Task 3. Derivation of preliminary wavefront c haracteristics 

Seidel or Zernike polynomial contributions were estimated after the polynomial- 
fit routine was integrated into the code. Contour maps of the pupil-plane phase 
were used to clarify the results. Errors in the wavefront reconstruction by phase 
retrieval were estimated. Preliminary conclusions were drawn on the conic 
constant error of the OTA primary mirror. 
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3.0 Technical Approach 

A direct (noniterative) technique for retrieving the phase distortion in the optics 
of instruments such as the Hubble telescope has been developed, and was 
used with Hubble telescope data. This technique took the digital information 
from distorted images of a star at two or more focus settings fairly close together; 
for example, at center focus and at ,say, half the depth of focus. The 
requirement was to obtain the local rate of change of the intensity at each pixel 
with respect to the amount of defocus, in each wavelength band. The more 
accurate this estimate was, despite noise in the image (including quantization 
noise due to limited numbers of bits per pixel), the more accurate was the phase 
map that would be generated. 

A simple partial differential equation (based on R. A. Gonsalves' work on 'phase 
diversity'), using this intensity data, was solved for the phase associated with 
the aberrated image. The complex field in the image plane was then Fourier 
transformed to obtain the complex field in a pupil plane. Fresnel wave optics 
propagation to any other pupil plane was done easily, when needed. The 
phase of this transform field was displayed then as a contour map to visualize 
the aberrations This phase was analyzed in terms of the usual polynomials like 
spherical aberration, astigmatism, and coma for convenient analytic 
characterization. 

This was the stage of work (Phase A) that identified the nature of the aberration 
for hardware fixes to the instruments (JPL’s WF/PC, for example) in order to 
compensate for the phase error in the telescope optics. Note that amplitude 
correction as well as this phase correction also could be applied at the pupil 
plane. Phase retrieval in each wavelength band should be consistent with the 
same optical path difference owing to the surface figure error of primary and/or 
secondary mirrors. With the proper variety of PSF's obtained with the WF/PC 
and the Faint Object Camera, with differing angular offsets, the retrieved phase 
maps could reveal the analytic form of the surface figure for the primary and 
secondary mirror, a necessary first step in the optical redesign task for the 
cameras. 

The following is a description of the direct phase retrieval method that utilizes 
the photon energy transport equation. This equation is applied in the vicinity of 
the image plane (near paraxial focus) but the equation is applicable at any 
general position. 
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3.1 Direct Phas e Retrieval Method 

FOURIER RELATIONSHIP BETWEEN PUPIL AND IMAGE PLANES 

Consider the simple geometry of a pupil plane on one side and an image plane 
on the other side of the lens equivalent of the telescope and camera system. F 
is the effective focal length. The wavevector magnitude is 


k- — 
k " A 


(3.1) 


The image plane field h(x.y) is obtained from an integral solution of the Fresnel 
wave equation, when the pupil plane field H(u,v) is given. 

h(x,y) = -iAFexp[ik(x 2 +y 2 )/2F)] 

| duj dvexp[-2/riB(u 2 + v 2 )]exp[-2/ri(ux + vy)]H(u,v) (3 2) 


where the spatial frequency components u, v are 
u = X/AF v = Y/AF 

and X, Y are the coordinates in the pupil plane. Because of the integration over 
u and v instead of X and Y, the coefficient in the expression for h(x,y) is -iAF 
instead of the well-known -i/(AF). The defocusing parameter B is proportional to 
the despace from the paraxial focus position. 

Omit the first phase factor (and coefficient) according to the R. Gonsalves' 
approximation, 


h(x,y) = J duj dvexp[-27riB(u 2 + v 2 )]exp[-2;ri(ux + vy)]H(u,v) 
This field obeys the following Fresnel diffraction equation 
-2ffl^ = V 2 h 

<9B 


(3.3) 


(3.4) 


where 


V 2 = 


+ 


<9x <9y‘ 


and where h = |h|exp(i</>) 


(3.5) 
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Eq. (3.4) is still valid even without the above Gonsalves approximation, when 
the distance to the focal plane F is replaced by the variable (z+F), and the 
defocus parameter B is replaced by BF/(z+F) in Eq. (3.2). The approximation 
holds for small z compared to F. 

The image field equation, which is complex-valued, yields two real-valued 
image equations, one for the derivative of the intensity with respect to B, and the 
other for the derivative of the phase with respect to B. 

These results apply to any plane, not just the image plane, except where a 
Fourier transform relationship to the pupil plane is invoked. 

FIRST IMAGE PLANE EQUATION 

The first equation, for conservation of energy, is (generalized Gonsalves case; 

R A.Gonsalves, "Phase retrieval by differential intensity measurements," J. Opt. 
Soa Am. A, 4, 166-170, 1987; also see K. Ichikawa, A. W. Lohmann, and M. 
Takeda, "Phase retrieval based on the irradiance transport equation and the 
Fourier transform method: experiments," Applied Optics, 27, 3433-6, 1988) 

(3.6, 

where the right-hand side parentheses contain the energy flux. All gradient, 
divergence or curl operations are with respect to the variables x and y, and not 
z. This equation holds generally in any plane, within the Fresnel propagation 
approximation, when B is replaced by its equivalent that is proportional to the z- 
variable. 

This energy transport equation is obtained from the field equation by separating 
the real and imaginary parts of the equation after substituting for the field its 
equivalent magnitude times phase factor. The transport equation comes from 
the imaginary part, while the second equation, given later, comes from the real 

part. 

The change AB in the defocus parameter B is proportional to the amount of 
defocus Az. 

AB = -Az ( 3 . 7 ) 

A value of defocus distance 

F 2 

Az = 8A^ 
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yields one wave change at the pupil edge (diameter D). The effective F no of 
the telescope and camera system is F/D. This particular length Az is the spatial 
period for the spots of Arago on the optical axis. 


The general form for the image-plane phase, <p , is obtained from the flux 

|h| 2 V0 = Vy/ + VxA (3.8) 

in terms of unknown potentials t/zand A. Note that d/V(curl A) = 0. So that the 
curl part of the flux does not help balance the change in intensity along the z- 
direction. 

The most important step is to solve Poisson's equation for V 


-K 


at 

oB 



(3.9) 


Note that V 2 y/ relates to the intensity change axially alone, while V 2 0 relates 
to the intensity change along a ray, axially, and sideways as well, as given by 
the direction of Vy/. 


Let it be noted that without the curl A contribution the contours of constant 
phase are parallel to the contours of constant intensity. 

SECOND IMAGE PLANE EQUATION 

The second equation that obtains from the field equation 


_27ri^- = V 2 h 

dB 


(3.10) 


is 


d(p 1 

Tt^r = - 

dB 2 




(3.11) 


the equivalent of the Eikonal equation. 

A spatially-varying refractive index term is included customarily. 

Because the unknown components of the gradient of the phase, in both axial 
and transverse directions, appear in this equation, there is no application for 
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this Eikonal equation as such. Of course, the original, simpler Fresnel equation 
for the complex field would be valid at any point and would be preferable to 
using Eq. (3.11) anyway. 

PHASE RETRIEVAL EQUATION SOLUTION 

The left-hand side of the phase retrieval equation is approximated by means of 
the small differences between two images that are slightly defocused from each 
other. 


= (Pixel-by-Pixel Differences) / AB 


* 


The defocus distance, that is, the difference of position A z is related to the 
difference in the so-called Goddard position as 

Az = 110 (1 .25) 2 x Goddard position difference. 

For example, Goddard position difference = 5 pm yields Az = 0.86 mm. 

SPECIAL IMAGE PLANE CYCLIC FLOW SOLUTIONS TO ADD 

Generally, the contours of constant phase are not parallel to the contours of 
constant intensity. Aside from trivial differences owing to an overall tilt to the 
phase wavefront, this difference of contour maps can be attributed to the curl A 
contribution, referred to here as cyclic flow solutions. This is in analogy to the 
vortex flows in hydrodynamics. Since the ±2 k phase jumps caused by these 
wavefront dislocations, do not affect the actual values of the fields themselves, 
the propagation to the near field or the pupil plane is unchanged as well. So 
the inclusion of these cyclic solutions in the image planes is not necessary for 
retrieving the phase in the pupil plane. However, these dislocations or vortices 
occur wherever the intensity is zero. This is especially true of the exit pupil 
plane due to the effect of the obscurations. 

The following discussion of the image plane is to better understand these cyclic 
solutions in general. The same equations for the field amplitude and phase 
hold for any plane, not just an image plane. And in particular, it is these same 
cyclic flows, which occur in the pupil plane, that are specially treated to yield the 
underlying smooth wavefront to be approximated by Zernike polynomials, for 
example. 

To represent straight vortex filaments parallel to the optical axis, the vector 
potential A is allowed to have a z-component only, as an approximation. The 
gradient operator applies here only to the x and y coordinates. Then, 

A = a(x,y)e 2 
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and 


Vx A = Vaxe z (3.12) 

and the added phase <po is given by 

|h| 2 V 0 O = Vaxe z (3.13) 


4^e z = VxV0 o = 



+ e z 



where the vorticity, in a discrete representation, is 

f = 2q i 5< 2 >(r-r i ) 

j 

upon allowing negligible variation in the z direction. This is an allowance that 
should be reexamined in the light of the focal plane in the classic case of the 
circular aperture discussed later. But for now, one obtains 



except at points n where the intensity |h| 2 vanishes. At these nulls the right- 
hand side of Eq. (3.14) is proportional to a delta function such that integrating 
over it yields a ±27i change (qj = ±1 )■ Let 



where, because of Eq. (3.14), 

V 2 * = 0 (3.15) 

The contours of constant a or x are parallel to those of constant intensity or y/. 
Laplace's equation is solved with the boundary condition that the circulation 
singularity strength be +2 k or -2ti wherever the intensity has a simple zero. 
These singularities are wavefront dislocations (like vortex filaments or sheets). 
With an explicit z dependence reintroduced, these vortex filaments and sheets 
may close upon themselves in a three-dimensional manner. This matter is 
pursued further in the section on special solutions. 
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Thus, 

X = ln|h| + ... 

where the additional non-singular term (denoted by the ellipses) is another 
function of \h\ that allows x to satisfy Laplace's Eq. (3.15) with the appropriate 
geometrical boundary conditions, and to follow the contours of constant 
intensity. 

The phase contribution <P o is given by 


v 2 0 0 = 0 


(3.16) 


where <t > o is the imaginary part of the analytic function (in a region avoiding the 
dislocations and made singly-connected by branch cuts), whose real part is X . 

These considerations of phase dislocations and vortex filaments apply to any 
plane, not just the image plane. In particular, the pupil plane phase has several 
dislocations, which are locally averaged out (in the Roddiers' POLE routine 
described later) to yield the desired smooth Zernike polynomial representation 
of the wavefront error and, hence, the primary mirror deformation. 

APPROXIMATE SOLUTION FOR IMAGE PLANE PHASE 


Solve for the phase from 


V0 = 


^ + V*xe z 

M 


(3.17) 


Note the mathematical identity that 

V • Vx x e 2 = 0 


Then, the phase is 


0 = V" 2 V-( 



where 

Vy = -/rV -2 V 


oB 


(3.18) 


(3.19) 
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The operations of inverse Lapiacian V 2 and so on are easily implemented in 
Fourier space, or in the case of cylindrical symmetry, in appropriate radial- 
coordinate integration. 

When (p 0 is unimportant, phase contours are parallel to intensity contours. 
Phase contours are sensitive to details of log (intensity). 


The contours of constant o perpendicular to the contours of constant 
intensity, in this approximation. The intensity contours are like streamlines to 
the flow, which is given by the gradient of this added phase. This is cyclic flow 
around the zero intensity points and confined by the intensity contours. 

With the phase determined, perform the Fourier transform of 


|h|exp(i</>) 

Obtain 

H = |H|exp(iO) (3.20) 

Then, obtain the exit pupil phase Basically, this is the phase to be retrieved. 

When the corresponding phase is obtained for the unaberrated case, that 
phase can be subtracted from this exit pupil phase to obtain the contribution 
due to the aberration itself. The amplitudes can be corrected also. 

SPECIAL VORTEX SOLUTION 

The classic case of the field distribution due to an illuminated circular aperture 
shows that the z-dependence of the A potential can be abrupt, contrary to the 
assumptions so far. The vortex filaments consist of a set of infinitesimally thin 
concentric rings, which lie in the focal plane alone and which encircle the 
optical axis. These lie at the locations of the Airy dark rings. Meanwhile, on 
axis are infinitesimal vortex rings, similar to dipoles or doublets, that lie at the 
periodic minima, the dark spots of Arago. Between these are the more famous 
bright spots. All these vortices, i.e., the rings and the doublets, are special 
because they have no extension in the z-direction. 

With a finite frequency spread and with the axial uniform flow superposed, some 
of the streamlines of energy flux circulate around these vortex rings if they are 
very nearby to any part of a vortex ring. But much of this flow goes by without 
making the trip around. There is a surface which divides the circulating flow 
from the noncirculating flow for each vortex ring. There are stagnation points for 
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the flow on these dividing surfaces. These points are intensity maxima, the Airy 
rings and Arago spots. See Fig. 3.1 for a sketch of this vortex ring configuration 
for the focal region for a uniformly illuminated circular aperture. 

The rings have circulation ail in the same sense of rotation around the axis, 
while the doublets have their rotation all together in the opposite sense. This 
sums to zero total angular momentum. 



Fig. 3.1 . Schematic of the streamlines of photon flow around the focal plane 
vortex rings (only two are shown) corresponding to the Airy dark rings and the 
whirls associated with one of the on-axis rings, corresponding to a dark spots of 
Arago. The stagnation points in the flow are denoted by s, where the bright 
spots or bright rings are located. The whirls are actually distributed uniformly 
around the circumference of each of their underlying rings. The geometrical 
focus is at the center of the Airy rings. Only parts of the nearly horizontal 
streamlines are shown. 

The width of the circulatory regions that surround each vortex ring depends on 
the frequency bandwidth of the light, according to J. F. Nye and M. V. Berry, 
Proc. R. Soc. Lond. A. 336, 165-190 (1974). Their Fig. 10 shows an enlarged 
drawing of the region in the immediate vicinity of a vortex filament. A section of 
a vortex ring would resemble a vortex filament. Their Fig. 10 is oriented at 90 
degrees to Fig. 3.1 . The intensity nulls, where the vorticity lies, can be thought 
of as caused by interference between the wavelets that diffract from the edges 
of the circular aperture, and by the interference at the focal plane between the 
shrinking cylindrical wave in front of the focal plane and the expanding 
cylindrical wave behind the focal plane. 

TRW 11 


In Principles of Optics, 6th ed.,Pergamon Press, Oxford, 1980, pp. 445-6, M. 

Born and E Wolf describe the phase distributions in the focal region. In their 
Figs 8 45 and 8.46, the little open circles denote these nulls. These are also 
the cross-sections of the vortex rings in the focal plane. Fig. 8.45 also shows 
the on-axis ring, which is infinitesimal in diameter and equivalent to a doublet, 
as in fluid dynamics. The ray paths shown above in Fig. 3.1, run perpendicular 
to the co-phasal surfaces plotted by Born and Wolf in their two Figs. 

Baranova et al. (Sov. Phys. JETP ,53(5), 925-929, May 1981 and J. Opt. Soc. 
Am. 73(5), 525-528, May 1983) discuss the vortex filaments, i.e., the wavefront 
dislocations that extend in the direction of the optical axis (z axis). Noisy 
wavefields cause speckle patterns and cause creation and annihilation of 
filamentary dislocations in pairs of opposite sign. The intensity is zero on these 
filaments. 

The treatment of the vortex filaments in Eqs. (3.12-3.17) could be generalized 
by allowing the gredient, divergence , curl and Laplacian operators to include 
the partial derivative with respect to z, too. The original energy transport 
equation would include the contribution from the divergence of the flux 
component in the z direction as well. The three-dimensional Laplace's 
equation for phase p 0 would be solved, except for the singularities where the 
intensity is zero. The vector potential A would lie nearly parallel to the vorticity, 
which exists only where the intensity is zero, i.e, in a direction perpendicular to 
both the gradient of the intensity and to the gradient of this phase. Earlier this 
vorticity was introduced as a two dimensional delta function with a coefficient 
such that this phase will change by ±2 k around a filament (whether straight or 
bent around in a ring). More generally, the filaments can be treated as three- 
dimensional configurations carrying this unit of strength. 

Then, the desired cyclic solution of Poisson's equation is obtained from 


V« 0 (p) = V X |j^^(p') 

Ip -p 


(3.21) 


(Biot-Savart law) for a general vorticity geometry. The position vector p (and p') 
is in curvilinear coordinates, two of which lie within surfaces of constant intensity 
and the third is perpendicular to them. 

This phase reduces to the solid angle subtended by a simple closed three- 
dimensional circuit of vorticity. The solid angle jumps by 4jt upon passing 
through the circuit. This causes an apparent discrepancy by a factor of two 
when comparing the two- and three- dimensional cases, i.e., linear angle 
around a length of filament vs. solid angle when passing through a closed 
circuit. 

When dealing with thin filaments along the z-axis (£-axis), as in the earlier 
discussion for the two-dimensional case, Eqs. (3.14)-(3.17), the analytic function 
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(3.22) 


X(P) 


= f^4 L -Iq l 5 (2) (p i '- p i i ) = S c, iH p x- p i 

J #v — n : i 


p -p 


where the vorticity is given a discrete representation at the prescribed points 
denoted by the two-dimensional pj_j. These points include all boundaries of the 
aperture and obscured regions, where the intensity is zero and where filament 
images might be used instead of boundary conditions on the vorticity. This is 
necessary to achieve parallel contours of intensity and of x. or orthogonal 
contours with respect to <p Q . The coordinates in Eqs. (3.21) and (3.22) are 
orthogonal to contours of constant intensity. 

Based on the approximation in Eq. (3.22), the corresponding analytic function 
for the cyclic phase is 


^o(ri) = Zqi ar 9(Pi-Pj.i) 

i 

A conformal transformation from p x to r x completes the specification of the cyclic 
solution. 

3.2 Imaoe Interpolation 

The width of the mathematical square aperture in the pupil plane exceeds the 
diameter of the primary mirror by a factor, n a X/Xq , which allows interpolation 
between the grid points that represent the pixels when a Fourier transform to the 
image plane is performed. This factor corresponds to increasing the resolution 
by exactly n a , which equals two or four, for example, independent of the filter 
wavelength. The reference wavelength, Xo , is 0.508 pm (15.24pm / F n0 ) for 
the OTA/PC combination, whose effective F n0 = 30. 

The intensity at each pixel is now shared by the n a x n a pixels. These are 
weighted according to relative intensities in a model calculation for images from 
the aberrated OTA. The amount of aberration used here is given by the 
equivalent of a conic-constant error that is typical of present knowledge about 
the error. The distance of the image planes from the paraxial focus position is 
an important parameter for the model calculation. 

The image field is cropped, usually by a factor of two reduction. This is to avoid 
aliasing when Fourier transforming back to the pupil planes. Aliasing, which is 
due to undersampling, is the undesirable overlay of data from the apparently 
equivalent periodic regions that arise from the assumption of sine and cosine 
Fourier representations of the data. 

3.3 Phase Unwrapping 

The pupil phase is obtained from the arctan(lmH/ReH), where the ar 9ument 
the arctan function is the ratio of the imaginary and real parts of the field H. This 
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phase lies between -rc and k, by definition of arctangent. The continuous phase 
distribution is obtained from this generally discontinuous distribution by the 
subroutine POLE that was obtained from the Roddiers (as reported at the Image 
Inversion Workshop). 

The method of unwrapping the phase is to add or subtract 2k from the original 
phase distribution at appropriate places. An elementary loop is formed by the 
four adjacent pixels of each grid square of the pupil plane. Each phase 
difference that corresponds to a side of each loop is checked to see if it exceeds 
in magnitude the value k. If so, then the appropriate 2k is subtracted or added 
to reduce each such phase difference to less than n in magnitude. The new 
phase differences are then summed around each such loop. Each sum, in turn, 
is checked for zero value. If not zero then the center of the offending loop, 
called a pole, is assigned a +2 k or a -2k value, as appropriate to the sense of 
circulation. Note that discontinuities of phase by the amount k can still exist. 
These corresponds to a phase change by going half way around a vortex 
filament. These filaments are taken to be of unit strength, i.e., ±2k, and not 
integer multiples because such multiple filaments are believed to be unstable to 
splitting into unit-strength filaments (B. Ya. Zel'dovich, et al. in Principles of 
Phase Conjugation , Springer-Verlag, Berlin, 1985, pp. 79-84). 

These poles are joined by straight lines (one line per pole) into closest pairs of 
opposites in sign of circulation. Excess poles are joined to the nearest 
boundary points of the aperture. The line joining such a pair of poles is called a 
dislocation by the Roddiers. This is the branch cut for the multivalued phase 
function. The vortex filaments themselves might better be called dislocations, 
like the screw dislocations of solids. 

The poles themselves are like vortices for the vector field formed by the phase 
differences. As many as hundreds of poles are found in the retrieved phase 
distribution, without using the averaging over azimuthal angle. With the 
averaging, the number of poles drops to tens, for example. These poles may 
have extension in the axial direction to form vortex filaments. Then, the Roddier 
branch cuts would form sheets. 

The phase distribution is obtained by summing the original phase differences 
along horizontal rows (working from the bottom upward) until a dislocation is 
encountered. Then, the appropriate jump of ±2k is applied in the vertical 
direction across the dislocation. At the end of this phase unwrapping the pixels 
on the dislocations themselves are given the average values of their nearest 
non-dislocation pixels. This last step removes any vestiges of curl contributions 
in the pupil plane phase map. 

3.4 Fnd-to-End Wave Optics Model 

The exit pupil phase contains the effects of the optical train even when free of 
aberrations. To determine the effect of the aberrations alone, the ideal phase 
should be subtracted from the retrieved phase, and the ideal, pupil amplitude 
divided into the retrieved amplitude. Furthermore, the defocusing operation on 
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HST is done by displacement of the secondary mirror, which displacement 
introduces additional aberrations besides defocusing that are modeled as well. 

A calculation of the ideal phase and amplitude properties was done by Fresnel 
propagation along the OTA and PC axes. This follows the Code V prescription 
provided by JPL. An efficient Fourier transform technique is used to implement 
the propagation from element to element in the train, in the presence of large 
changes in magnification. This was related to the Talanov transformation, 
according to S. Ebstein, and used by the Roddiers in their repeated Fresnel 
transform (double Fourier transform) technique, as described at the first HARP 
Image Inversion Workshop. 

The mirror surfaces for OTA and PC are modeled by appropriate conics of 
revolution (hyperboloids). The field flattener lens is modeled as a simple 
paraboloid (quadratic phase factor) lens and spacer. 

Test runs have shown the image at the focus of the OTA, namely, at the 
pyramid, to be about as expected. However, the image produced at the image 
plane of the PC was not like observed images. The propagation algorithm 
through the PC was very similar to that for the OTA. The nature of the difficulty 
was not clear. 

While the end-to-end model was being debugged, phase retrieval test and 
interpolation results were obtained using a primary and a camera mask for 
specifying the input of relative intensities, i.e.,the zeros and ones of the two 
pupil functions. For the simple simulations and interpolations, these two masks 
were merely multiplied together and input at the primary mirror. However, the 
end-to-end model has the PC6 mask input at the relay (camera) primary (or 
somewhat upbeam), while the OTA mask was input at the OTA primary. These 
masks consist of the near-central obscurations, the spider-support obscurations, 
and, for the OTA, the three low-reflecting mounting pads. The PC6-channel 
orientation of these elements was considered only. 

3.5 Zernike Polynomial Fitting 

The pupil plane wavefront phase map was expressed as a sum of polynomials, 
chosen to be the orthogonal set of Zernike type with a 0.33 fractional radius for 
a central obscuration. The first 22 of these were tabulated in the HST OTA 
Handbook, vers. 1 .0, May 1990. These were normalized to unity integral over 
the clear aperture area, and, hence, unit rms value. 

The least-squares routine obtained the coefficients of the polynomial expansion 
by fitting the phase factor, the exponential function of i times the phase, rather 
than the phase function itself. In part, this circumvents the ambiguities of 
multiples of full waves. However, phase unwrapping was done to remove these 
ambiguities anyway. 

The error metric that is to be minimized is 
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(3.23) 


Error = £|exp(i<I> k )-exp(i@ k )| 2 |H k | 2 / £|H k 

where O k is the digital data to be fit and 0 k is the sum of Zernike polynomials. 
The subscript k denotes data at points x|< within the annular aperture. The 
quantity H k is the retrieved pupil field distribution. The sum with respect to the 
subscript k in the denominator is the total power within the annular aperture. 

@k=E a M Z M ( x k) (3.24) 

The normal equations for minimization are used to solve for the increments Aa M 
in the polynomial coefficients. 


i^ = -2£sin(4> k -0 k )Z u (x k ) 

4^U2£cos(<t> k -0 k )Z M (x k )Z N (x k ) 

_1 ^Error 

MN ^ a N 


Aa M = 


<9 2 Error 

dada 


(3.25) 


The matrix element is the MN th element of the inverse matrix of second 
derivatives of the error metric. Sums are over the repeated index on the various 
quantities on any of the right-hand sides. 

These equations are nonlinear. During an iteration procedure they depend on 
the previous values of the polynomial coefficients. Schemes for improving 
convergence of the iterations that are performed in this minimization routine are 
straightforward but untried here. Without the smoothing effect of phase 
unwrapping this simple iteration procedure is usually unstable. 

3.6 Error Estimates for Phase Retrieval 

Noise Effects 

By means of a reverse Strehl ratio calculation, i.e., averaging over the image 
plane instead of the pupil plane, and by the Cumulant Theorem (see M. 
Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, U.S. Govt, 
Wash., DC, 1972, p. 928 (item 26.1.12); and R. Kubo, M A Stochastic Theory of 
Line Shape and Relaxation," in Fluctuation, Relaxation and Resonance in 
Magnetic Systems, ed. D. Ter Haar, Plenum, NY, 1962, pp. 23-68) on 
exponential functions, the variance of the pupil plane phase (in radians) is 
obtained from 
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(exp(i<50)) = exp(i(<50)--(<S0 2 )) 

= (exp(5ln|h| + i<5</>)) = exp(i(<5ln|h|50) + -[((5ln|h|) 2 ^)-(<5</> 2 ^]) ^ 2 g) 


where the usual factor exp(ln|h|) was included in the definition of the averaging 

in order to create the aperture function by Fourier transform, and was not 
included as part of the contribution to the phase. The imaginary-valued cross- 
correlation term contributes to the average pupil phase. 

Then, 

(54> 2 ) = -^((Sln|h| 2 ) 2 ) + (50 2 } (3.27) 


( ) denotes the spatial weighted-average over image area, instead of over the 
usual pupil area for the usual Strehl ratio. The 8 denotes the statistical 
deviation due to noise or other random effects. Specifically, 

f d 2 rexp(-2/riu • r)h 0 (r)exp(<5ln|h| + i <50 ) 

(exp(i<50)) - Jd 2 rexp(-27riu-r)h 0 (r) 


where 


J d 2 rexp(-2;riu • r)h o (r)f(5ln|h|,<50) 
Jd 2 rexp(-2;riu-r)h 0 (r) 

for any function f of the fluctuations. The underlying image field h 0 produces the 
aperture function by Fourier transform. 

It is simple and reasonable to assume that the statistical probability distribution 
is Log-Normal for image intensity and Normal, i.e., Gaussian, for image phase. 

The subtractive term on the right-hand side of Eq. (3.27) will be smaller than the 
second term on the right-hand side, according to the estimates below, and as 
required by the positivity of the left-hand side of that equation. 

The fluctuation in logarithm is related to the small noise fluctuation by 
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<5ln|h| 2 = ln(signal + noise) - (ln(signal + noise)) 

= <5noise / (signal + (noise)) 

so that, for example, 

((Sln|h| ) 2 ) = (1 + SNR) 2 (3.2! 

for a simple exponential noise probability distribution over the image area, for 
which 

(<5noise 2 ) = (noise) 2 


SNR denotes the signal-to-noise ratio for area-averaged quantities. This does 
not refer to temporal statistics at a given pixel, which are close to Poisson 
statistics, for which 


(<5noise 2 ) = (noise) 

The Cumulant Theorem was applied , over the image area, to assumed spatial 
Gaussian statistics for the random variables in the exponents. This simple case, 
for which the cumulant sum is truncated at the second term, is the limit for any 
type of statistics when the fluctuations are weak and/or compact. That is, the 
value of the product of rms modulation of spatial gradients of the random 
variables and correlation distance is smaller than one. If the correlation 
distance is only on the order of a pixel or two, then the rms gradient of 
In(intensity) need be small over that distance. 

But with the more relevant Poisson noise statistics over the image area, in 
terms of numbers of noise equivalent-photoelectrons and signal photoelectrons 
per image, Eq. (3.29) becomes 


/ , l2 o\ (noise) (noise) 1 

( (5lnlh| ' ) = (signal +<noise)) 2 ~(1 + SNR) 2 


(3.30) 


SNR of greater than tens are found for the cropped images, while at least a few 
hundred noise equivalent-photoelectrons per pixel would be found. About 1 28 
x 128 pixels are used per image. 

The image phase, <p, excludes tilt and focus terms, in effect, because of the 
Fourier transform phase term: -2jtu-r, which is implicitly included. The quantity 
u is the pupil plane position, in spatial frequency or in spatial pupil-plane 
position divided by \F, at which the variance of the phase <D is being estimated. 
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With no central obscuration, this would be at the center of the pupil plane 
corresponding to u=0, in the spirit of a Strehl ratio. However, the variance in 
Eq. (3.27) is evaluated at a typical point in the clear aperture, where the 
aperture function is unity. 

An estimate of the phase error is obtained by way of the phase retrieval 
equation, Eq. (3.18), when a simplification is made that the term involving the 
gradient of the reciprocal intensity can be neglected. This is explained as 
follows: 

In the absence of dislocations, near an intensity maximum, 


<f> = V" 2 V-( 



) 


where 
Note that 

y2 =1 j. r — + — 

r dr dr r 2 dd 2 

for the two-dimensional case dealt with here. For angle-independent quantities, 



where the lower limits of integration were chosen to eliminate logarithm and 
constant terms. 

Upon making use of the operator identity 

V-V~ 2 V = 1 

which follows from interchanging the order in which partial derivatives are done, 

V~ 2 V = VV~ 2 

and from using the definition of the Laplacian operator 

V-V = V 2 
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and, for the sake of simplicity, upon neglecting the gradient of the reciprocal 
intensity near an intensity maximum, where the gradient vanishes, the retrieved 
phase is approximately 



-Lvv- 2 v 

|hf 


<?B 


|h| 2 <?B 


or 



<9ln|h| 2 

<9B 


( 3 . 32 ) 


Then, with the finite difference approximation to the derivative with respect to 
the amount of defocus between the two image planes 

<?ln|h| 2 __ Aln|h| 2 
dB " AB 

and with the variance of this quantity 



<?ln|h| 

<?B 


2 Yl 


J) 


((Mln|h| 2 ) 2 )_2((5ln|h| 2 ) 2 ) 
(AB) 2 = (AB) 2 


and, with the variance of ^ 2 MN 
given by 


J 2 k s+n 




itsaiiq 

J 2 K J 2 K \ s + n s + n , 
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then 


jal {(gn) 2 )dg lx .i ^odp i>9 i 

= 2(2 re) 2 (s + n) 2 “ 2(2jr) 2 (1 + SNR) 2 (noise) 




<9ln|h| 

~^B" 


2 A 



af 2 0 1 

X 2 AZ 1 + SNR 


n 


pixel 


(3.33) 


where A z is the defocus distance, and 7tro 2 = No 2 dpjxei 2 .where ro is the radius of 
the image region. N 0 2 is the number of pixels utilized per image. There is an 
extra factor of two in the numerator because the variance is of the difference of 
two fluctuating ln|h| 2 quantities in Aln|h| 2 . The 5 denotes deviations from the 
appropriate spatial mean. The pixel size dpj X6 | = A<)F no Note that the subtractive 
term in Eq. (3.27) can be ignored when relating Eq. (3.33) to the pupil phase 
variance. 


Slipping the angular brackets inside the integrals is done by an appropriate 
interchange of the order of integrations. The statistical approximation for the 
noise correlation function 


{5n(r’)5n(r")> = {(5n(r')) 2 )d 2 xel 5 (2) (r'-r") 


(3.34) 


was used. This is equivalent to treating the noise in the different pixels as 
statistically independent of each other. This whole procedure is somewhat 
cavalier, but useful for a rough estimate of the variance. 

The average signal density per unit area was denoted by s and the average 
noise areal density by n. That is, 


n = (noise) /(/rr 0 2 ) 

Note that the noise per pixel is 

"pixel = ( n0ise > 1 N 0 
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This phase error is inversely proportional to ^ )^jp^xa\ with SNR 

between 10 and 100, with X =0.959 X 0 at 487 nm, with N 0 = 128, and with Az = 2 
A F no 2 , as is the case for these data, the rms phase error is between 0.01 and 
0.001 waves, respectively, when there are a typical 100 equivalent noise 
photoelectrons per pixel, a convenient reference level. 

Defocus Error 

Other sources of error might dominate this noise error. The so-called breathing 
in and out of the focal position is such a case, especially if the characteristic 
time is comparable to the data integration time, while the amplitude of the effect 
is comparable to the separation distance of the pair of image planes used for 
phase retrieval. There is no indication for such a difficulty from these results. 

Shorter exposure images, with the narrower central peaks, were dealt with first. 
These would be less affected by the breathing defocus action and also less by 
alignment jitter. 

If the error in the image plane separation Az is important, then 



Here the angular brackets denote averaging over the ensemble of fluctuations 
in the defocus distance. 

For simplicity in evaluating the second derivative, the image intensity is 
temporarily assumed Gaussian-shaped, having width that fits to the diffraction 
from a uniformly-filled aperture and that has blurring due to defocusing a 
distance z from paraxial focus. The blurring due to spherical aberration is 
deferred from consideration for the moment. 


In|h | 2 = ln|h 0 | 2 


/rr 


l2AF no J 
1 + w 2 


and 
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|2 ..,2 o .«.2 n ( •— A 


crln|h| 0 w 3w -1 
<9B 2 =2 B f (1 + w 2 ) 3 


/rr 


2AF^ , 


where 


W = 


7rB 


7TZ 


8A 2 F 2 0 . 16AF, 


no. 


(3.36) 


Then 



f- 


2 ln|h| 2 tit 2 1 
<® 2 16 , 


^([5(AB)] 2 ) 


(3.37) 


Note that the rms error increases with the fourth power of the radius in the 
image, when the quadratic dependence in the second derivative is included. 


A fractional error 5(AB)/AB of 10 percent (say, 0.5 pm out of 5 pm in Goddard 
units) yields 8(Az) = 0.0086 cm. When z = -1 .96 cm (at -114 Goddard units) and 
the quantity for blurring w = -8.8, then an rms phase error of 0.04 waves is 
estimated at the edge of the image. If the variance is averaged over the image 
area, the final rms error relevant to the pupil phase is a factor 2.236 smaller. 
Then, the pupil phase error is approximately 0.018 waves. 


These are only rough estimates, while errors in AB are not yet known. There is 
no strong evidence that the breathing motion of focus jitter has greatly affected 
the difference quantity AB in the measurements at 487 nm that were used here. 


Other Errors 


The quantity B itself is to be measured from the paraxial focus position, which is 
a parameter that is varied as part of the least squares fitting to Zernike 
polynomials. 

Another error source has been the unwrapping of the phase, where the pupil 
wavefront dislocations are smoothed over, in a final step. The full two- 
dimensional Fourier method of phase retrieval yields a wavefront with many 
dislocations. Zernike polynomial fits to these smooth unwrapped phase 
distributions seemed consistently shifted from the corresponding fits to the 
azimuthally-averaged (one-dimensional) phase distributions. This shift has not 
been made quantitative yet, due to the limited number of runs in the two- 
dimensional case, and due to the wider polynomial set allowed for this case. 

The Zernike fitting introduces errors due to a dependence on the number of 
polynomials that are allowed in the set to be fit, even in the simple cylindrically 
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symmetric case. When the HST OTA Handbook Z 22 polynomial (fifth order 
spherical aberration) is included along with the usual Zi (piston), Z4 (focus) and 
Z11 (third order spherical aberration) , a limiter on the amount of change in 
coefficient for Z 22 at each iteration sometimes has to be introduced to prevent 
numerical instability. Apparently, there is some underlying interdependency of 
the variables that makes the inverse matrix elements small when involving this 
last polynomial. 

Error due to digitization with 12 bits has not been a difficulty. The dynamic 
range of 0 to 4095 seems adequate to handle the large central peak of the point 
spread function. Saturation effects on the images have been avoided by 
choosing those image pairs with the narrower central peaks. 

Shot noise and read-out noise were the main contributors to the noise in the 
signal-to-noise ratio estimated to be between 10 and 100. 

Flat-fielding error and errors due to cosmic-ray spike removal were not 
evaluated. These were not expected to contribute more than the shot noise in 
the image central region that contributed most significantly to this phase 
retrieval. 
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4.0 Major Findings 

4.1 Task Results 


4.1.1 Task 1 . HST Point Source Images 

Digital data of point source images were obtained from the JPL VAX facility via 
magnetic tape and SPAN in September. The data were of the same source at 
different wavelengths and at different focal positions. The focal position was 
changed aboard HST by movement of the OTA secondary mirror. Attention was 

mainly directed to the data set designated PC6F487N_G1 .FITS, ...HI II .... 

and ...G2 H2. ...... 12.... These data correspond to narrow band 487-nm signal 

with the OTA secondary mirror at the Goddard positions of -5, 0 and 5 |im. This 
early choice was made on the basis of proximity of the positions to the 'best' 
focus, and the small despace between images. 

Digital data of point source images were obtained from the JPL VAX facility via 
magnetic tape in November. The data were of the same source at different 
wavelengths and at different focal positions. Attention was directed to the data 
set designated PCF487N_01&2.FITS and PCF487N_P1&2.FITS. These data 
correspond to narrow band 487-nm signal with the OTA secondary mirror at the 
Goddard positions of -267 and -260 ^im. At these same Goddard positions, but 
with narrow band 889-nm signal, the set PCF889N_01 &P1.FITS was also 
analyzed. These choices were made on the basis of small despace between 
images, while representing cases of severe blurring. 

The images appear to require no additional registration within the image frames 
for centering the images on each other. 

4.1.2 Task 2 . Phase Retrieval Code 

The phase retrieval code was implemented on other simple test cases. Large 
arrays (256 X 256) were used to reduce the raggedness of the circular edges 
on the numerical grid. Effects of filtering out high spatial frequencies in the pupil 
phase function was explored. Flat phase to less than 0.06 waves was retrieved 
for the input flat phase case with the centrally-obscured circular aperture. 

The phase retrieval code was modified to deal with a cylindrically-symmetric 
average of the data. The differential equation of this direct phase retrieval 
technique is now integrated in the radial coordinate in the image plane. This 
substitutes for the Fourier transform equation-solving technique used for the 
general case. 

The rings of zero intensity in the image plane caused half-wave phase jumps. 
These are dealt with more accurately in the cylindrically-symmetric, non-Fourier 
equation-solving case. These zeros arose from the hard edges of the apertures 
in the near-field. The reconstruction of the pupil plane amplitude and phase 
would better show the obscuration characteristics of telescope and camera 
when the half-wave phase jumps are calculated better in the image plane 
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phase distribution. However, an accurate reproduction of the pupil plane 
characteristics would require the full two-dimensional phase retrieval technique. 


The energy transport equation for the cylindrically-symmetric case is 


_4 = 4 [r|( r, B) fl 

<2B r<?r dr 


(4.1) 


where l(r,B) is the image intensity distribution (averaged over azimuthal angle) 
in the plane, whose axial coordinate is proportional to the defocus parameter B. 

This equation is integrated to solve for the phase in the image plane. 



1 

r'l(r\B) 



<?l(r M t B) 

dB 


(4.2) 


This solution has the boundary conditions built in that the phase is zero at the 
center and that its slope is also zero there. 

Note the division by the intensity within the last integration. The zeros of 
intensity, as discussed before, lead to half-wave jumps in phase. In principle, 
these jumps would result from the integration. However, in practice the values 
of intensity in the vicinity of each zero are known only very roughly on the image 
plane grid. The intensity values are small near the zeros and introduce errors in 
the result. For reasons that involve reintroducing terms that were dropped in the 
paraxial or Fresnel approximation right at the beginning, this denominator was 
replaced by the following more accurate and less offending quantity: 


Va 2 +b 2 

where the original denominator 

a = r'l(r\B) 

and the new contribution 


Jo k<9B 


where k =2 k/X. 

The quantities a and b are proportional to the components of the local ray 
direction. This follows from a good approximation that the rays lie on the 
surface of constant enclosed power, which is the area integral of the intensity up 
to a radius r. It follows from the energy transport equation that while the 
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curvature of the rays is slight so that v2 0 = 0 as in incompressible flow, (the 
Laplacian now includes the second order z-derivative), and while the rays are 
still mainly parallel to the z-axis and the intensity changes are mainly radial, the 
slope of a ray is approximately 

dr _ dPJ_dz 

dz _ _ <?P / dr (4.3) 

where the enclosed power P is 

P = 2^drri(r',z) (4 4) 

Note that z and B are proportional. Note also, that exactly at a vortex position, 
the quantity b, as well as a, is zero, whereby the ray direction is indeterminate. 

However, with this modified denominator, the phase <p is correctly given by the 
ray path integral along the radial direction, with the appropriate sine of the 
angle of this ray with respect to the axis, instead of the tangent of this angle. 

The two trigonometric functions are usually nearly equal, except where the 
intensity is small, and where the dominance of the axial wave vector component 
is no longer a good approximation. Nevertheless, the new integrand is now 
bounded, since the sine, in magnitude, is less than or equal to one. 

A similar revision of the denominator in the two-dimensional phase retrieval 
algorithm was introduced. 

The wave optics propagation routine for dealing with the different planes 
containing obscurations, apertures and other optical elements of the OTA and 
the PC was developed. Details of the pads and spiders were introduced for the 
OTA and the PC models. The secondary mirror position can be changed to 
cause defocusing of the image. The aberrations introduced by this change can 
be modeled with this code, as well as providing test images for phase retrieval 
and for interpolation between the pixels of actual data. 

As an instructive exercise, the prescription, provided by JPL, for the HST OTA 
and the PC, was entered into Code V at TRW. Various plots were made to 
clarify the configuration. 

Appropriate geometrical specifications were introduced into the wave optics 
code to nearly complete the end-to-end propagation capability that this code 
would provide. The strategy to distinguish the phase aberration contributions of 
the OTA primary from those of the secondary mirror would depend on this code, 
when used for cases of different incident field angles for the star images. Also, 
besides changing the focal position, motion of the secondary mirror produces 
small aberrations to the wavefront. 
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4.1.3 Task 3 . Preliminary Wavefront C haracteristics 

The test cases indicated that the polynomial fitting was sensitive to even very 
small effects such as the rapidly varying intensity and phase conditions at or 
near the pupil edge, and the discrete square grid of input data that impressed a 
four-fold symmetry on the retrieved phase. 

The even-simpler test case of uniform phase (i.e., no spherical aberration) over 
a circular aperture, with a circular central obscuration (like the OTA case), 
revealed similar four-fold symmetry features, even for the 256 X 256 element 
pupil-grid case. 

The (second) iterative fitting scheme outlined by J. R. Fienup (HARP Image 
Inversion Workshop, Nov. 1990) was implemented to decrease sensitivity to 
phase unwrapping. However, the normal equations for minimization were used 
to supply the changes in coefficients for the next iterations. 

The Goddard position for paraxial focus was varied from case to case until a 
best fit was obtained according to the least residual error. The interpolation 
between pixels was redone for each of these cases. 

The phase unwrapping routine of C. and F. Roddier, discussed briefly at the 
November HARP Workshop, was implemented. The routine has run, where 
hundreds of wavefront jumps (dislocations) were unwrapped and smoothed for 
the Fourier, two-dimensional phase retrieval technique, while only tens of 
dislocations arose and were handled in the cylindrically-averaged version of 
the phase retrieval. 

4.2 Separating Primary and Secondary Aberration Contributions 

The strategy for separating the contributions from the three sources: the primary 
mirror, the secondary mirror and the camera, was still being formulated. The 
fact that the phase contributions would be weighted differently with changes of 
field angle would be critical. 

The geometrical optics limit would yield the resultant OTA phase as the sum of 
the phase contributions from primary and secondary mirrors along any given 
ray path between them, and then on to the image. Wave optics alone would 
make a very small, angle-dependent, correction to this sum. 

h(r) = J (0'p(K)WK-AK)f s ([K + AK-^]M) 


where K is the Fourier transform wavevector variable in the two-dimensional 
integral, k is the usual 2 tc/X, AK is the wavevector corresponding to the off-axis 
field angle, M is the secondary magnification (based on the secondary's 
illuminated diameter), and F is the effective focal length of OTA 
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f p (K) = F.T.[exp(iAO p )] 


is the coherent point spread function for the OTA primary, where AO p is the 
phase deviation from paraboloid shape of the primary. 

f ps (K) = Mexp[-iK 2 Md / 2k] 

is the wave propagator between primary and secondary (distance apart, d). 
f s (K) = F.T.[exp(iA4> s )] 

is the coherent PSF for the OTA secondary, where A<D S is the phase aberration 
of the secondary. F.T. denotes Fourier transform operation. 

In the geometrical optics limit: f ps (K) = M so that 
h(r) = Mj ^f p (K)f,([K + AK-^]M) 

\dK) r 

= F.T.[exp(iAO p + iA<D s )] ( 46 ) 


where A<D p + AO s is evaluated at points on the same ray that goes from 
primary to secondary mirror and onward. 

The small correction for wave optics is obtained by Taylor expanding the 
propagator f ps . 

Ah(r) = M f f 0 (K){-i|K - AK| 2 Md / 2k)f s ([K + AK - ~]M) 

J (2 kY f (4.7) 

This correction term is inversely proportional to the Fresnel number of the 
secondary mirror, owing to the fact that the typical spatial frequencies 
associated with K are only of the order of the reciprocal diameter of the primary 
The two magnification factors lead to a Fresnel number based on the diameter 
of the secondary. This Fresnel number is very large, so that the effect of wave 
optics, compared to geometric optics is small in this case. 

The important point is that the sum of phases, primary plus secondary, is 
different for different angles of input wavefront tilt, given by AK. If the 
contribution from the primary is nearly the same then the change in the sum is 
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mainly due to the secondary. Thus, data for different off-axis objects would yield 
the required information on the secondary mirror. 

This procedure was not attempted, however. Effects of the images lying near a 
corner of the CCD chip rather than near the center, and other effects, needed to 
be evaluated first. 

4.3 Conic-Constant Error 

The relationship between the Zernike spherical aberration coefficient Z u (pm) 
and AK' the error in conic constant is related to the fourth-order dependence on 
p, the fractional radial distance in the pupil plane. For the spherical aberration 
coefficient, the optical path difference (in meters) due to the mirror deformation 
is 

OPD = AK’ r p 4 /(8 R 3 ) = 0.5 Z-, 1 1 0' 6 1 6.896 p 4 (4.8) 

where r = 1 .2 m, the radius of the primary mirror, and R = (1 1 .04 m)/(1 .2 m), the 
ratio of the radius of curvature to the mirror radius. The factor of 0.5 converts 
wavefront phase to mirror deformation. The factor 16.896 provides the correct 
spherical-aberration polynomial normalization, with a 0.33 obscuration ratio, so 
that over the clear aperture, the area integral of the square of the polynomial is 
unity. Then, the coefficient is the rms value of spherical aberration over this 
clear aperture annulus. 

Values for at the wavelengths 0.487 and 0.889 pm were -0.27 and -0.28 
pm. The reason for the discrepancy is not known. These two wavelength cases 
were from HARP 1 A and 1 B data, respectively. A similar difference of Z-n has 
been noted by R. Lyons (at the second HARP workshop) for these data sets, 
irrespective of wavelength. 

At 487 nm, the Zernike polynomial coefficients (0.33 obscuration ratio) in waves, 
for piston, focus, third order and fifth order spherical aberration, with the paraxial 
focus distance at -1 14.0±0.02 Goddard units, were 

based on HI -II data: 

Zi =-3.69±0.01, Z 4 = -2.95±0.01, Z u = -0.56010.002, Z 22 = 0.0378±0.002; 
based on G1-H1 data: 

Zy = -3.5910.02, Z 4 = -2.8510.01 , Z u = -0.55310.003, Z 22 = 0.000610.003. 
Compare these with 

Zi = -3.54, Z 4 = -2.70, Z 11= -0.560, Z 22 = 0.0, 

for a pure quartic: 

-0.273 pm 16.896 p 4 

for the effect of conic constant error. Note that Z 22 should be negligible (of order 
10' 4 ) for the actual (hyperboloid) conic figure of revolution. The first factor in the 
coefficient of p 4 is the derived Z\ \ in waves times the wavelength in pm. 
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At 889 nm, the Zernike coefficients in waves, with paraxial focus distance at the 
later observation date approximately equal to -11911 Goddard units, were 

Based on 01 -PI data: 

Z , = -2.15±0.01, Z 4 = -1.6410.01, Zn = -0.319±0.005, Z 22 = -0.095110.0003. 
Compare with 

Zi = -2.01 , Z 4 = -1 .54, Zn= -0.319, Z 22 =0.0, 

for a pure quartic: 

- 0.284 pm 16.896 p 4 

for the effect of conic constant error. Again, the first factor in the coefficient of p 4 
is the derived Zn times the wavelength in pm. 
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5.0 Conclusions 

The retrieval results indicated sensitivity of the fitting of Zernike coefficients to 
symmetry anomalies in the data, such as the four-fold symmetry due to the 
square grid of pixel intensity data. The data, when averaged over angle in the 
image plane, and the direct integration technique, using only the radial 
coordinate, helped to overcome this difficulty. 

The phase unwrapping (full-wave jump smoothing) routine of the Roddiers 
handled the fragmented phase that was retrieved, especially when the input to 
the retrieval was image data that was averaged over angle around the central 
pixel. 

The derived Zernike coefficients depended on the paraxial focus position 
chosen. However, when the Goddard position of the paraxial focus was varied 
as a parameter, a definite best fit was found. This was characterized by a 
significantly lower residual error than found at neighboring values of paraxial 
focus positions. The paraxial focus was, thus, found to be -1 14.0 pm for the 
HARP 1 A 487-nm data. The results for the HARP IB 889-nm data were less 
definitive with respect to low residual error, owing to the much more blurred 
images and owing to some constraint needed on the Z 22 coefficient during 
fitting. 

The Zernike coefficient Z 11 (HST OTA Handbook) was -0.27 and -0.28 pm for 
the HI -II and 01 -PI pairs of images, at 487 and 889 nm, respectively. By 
using the formula obtained from the previous OPD relation, 

AK’ = Z 11 /22.802 (5.1) 

the error offset in conic constant is -0.01 18 and -0.0123. The rms error in the 
Zernike coefficient was ±0.005. This was estimated from the standard deviation 
for fitting (converged iterations) of the polynomial coefficients at 889 nm, 
associated with a narrow range of paraxial focus distances (-1 19±1 ). The error 
due to polynomial fitting in the 487 nm data was somewhat less. 

The design conic constant was -1 .0022985. Then, the estimated conic constant 
is -1 .01 41 for the 487-nm HARP 1 A data and -1 .01 46 for the 889-nm HARP 1 B 
data, with rms error of ±0.0002 due to polynomial fitting. 
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